"""
Grid Generation Utilities

Provides functions for creating uniform grids on L-shaped domains, including
coordinate generation, spacing calculations, and grid-related utilities.
"""

"""
    build_grid(n::Int, ::Val{2}) -> (Vector{Float64}, Vector{Float64}, Float64)

Build a 2D uniform grid on [-1,1]² with n interior intervals per axis.

# Arguments
- `n::Int`: Number of interior intervals per axis
- `::Val{2}`: Dimension specifier

# Returns
- `x::Vector{Float64}`: x-coordinates (including boundaries, length n+2)
- `y::Vector{Float64}`: y-coordinates (including boundaries, length n+2)  
- `h::Float64`: Grid spacing (uniform in both directions)

# Example
```julia
x, y, h = build_grid(32, Val(2))
@assert length(x) == length(y) == 34
@assert h ≈ 2.0/32
```
"""
function build_grid(n::Int, ::Val{2})
    x = range(-1.0, 1.0, length=n+2) |> collect
    y = range(-1.0, 1.0, length=n+2) |> collect
    h = x[2] - x[1]  # Uniform spacing
    
    return x, y, h
end

"""
    build_grid(n::Int, ::Val{3}) -> (Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64)

Build a 3D uniform grid on [-1,1]³ with n interior intervals per axis.

# Arguments  
- `n::Int`: Number of interior intervals per axis
- `::Val{3}`: Dimension specifier

# Returns
- `x::Vector{Float64}`: x-coordinates (length n+2)
- `y::Vector{Float64}`: y-coordinates (length n+2)
- `z::Vector{Float64}`: z-coordinates (length n+2)
- `h::Float64`: Grid spacing (uniform in all directions)
"""
function build_grid(n::Int, ::Val{3})
    x = range(-1.0, 1.0, length=n+2) |> collect
    y = range(-1.0, 1.0, length=n+2) |> collect
    z = range(-1.0, 1.0, length=n+2) |> collect
    h = x[2] - x[1]
    
    return x, y, z, h
end

"""
    grid_coordinates_2d(x::Vector{Float64}, y::Vector{Float64}) -> (Matrix{Float64}, Matrix{Float64})

Create 2D coordinate meshgrids for plotting and analysis.

# Arguments
- `x::Vector{Float64}`: x-coordinate vector
- `y::Vector{Float64}`: y-coordinate vector

# Returns  
- `X::Matrix{Float64}`: x-coordinates meshgrid
- `Y::Matrix{Float64}`: y-coordinates meshgrid

# Example
```julia
x, y, h = build_grid(10, Val(2))
X, Y = grid_coordinates_2d(x, y)
@assert size(X) == size(Y) == (12, 12)
```
"""
function grid_coordinates_2d(x::Vector{Float64}, y::Vector{Float64})
    nx, ny = length(x), length(y)
    X = repeat(x, 1, ny)
    Y = repeat(y', nx, 1)  
    return X, Y
end

"""
    grid_coordinates_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}) 
    -> (Array{Float64,3}, Array{Float64,3}, Array{Float64,3})

Create 3D coordinate meshgrids.
"""
function grid_coordinates_3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64})
    nx, ny, nz = length(x), length(y), length(z)
    
    X = zeros(nx, ny, nz)
    Y = zeros(nx, ny, nz) 
    Z = zeros(nx, ny, nz)
    
    @inbounds for k in 1:nz, j in 1:ny, i in 1:nx
        X[i, j, k] = x[i]
        Y[i, j, k] = y[j]
        Z[i, j, k] = z[k]
    end
    
    return X, Y, Z
end

"""
    compute_grid_metrics(h::Float64, Re::Float64) -> NamedTuple

Compute important grid-related metrics for stability analysis.

# Arguments
- `h::Float64`: Grid spacing
- `Re::Float64`: Reynolds number

# Returns
Named tuple with:
- `h`: Grid spacing
- `cfl_limit`: CFL stability limit (advection)
- `diffusion_limit`: Diffusion stability limit  
- `peclet_number`: Grid Péclet number
- `recommended_dt`: Recommended time step

# Notes
For explicit time stepping of vorticity transport:
- CFL condition: Δt ≤ h/|u|_max  
- Diffusion condition: Δt ≤ h²Re/4 (conservative)
"""
function compute_grid_metrics(h::Float64, Re::Float64)
    # Conservative factors for stability
    cfl_factor = 0.45      # Conservative CFL number
    diffusion_factor = 0.14 # Conservative diffusion stability factor
    
    # Stability limits (assuming |u|_max ≈ 1 for lid-driven cavity)
    cfl_limit = cfl_factor * h
    diffusion_limit = diffusion_factor * h^2 * Re
    
    # Grid Péclet number (measure of advection vs diffusion)
    # Pe = |u|h*Re (assuming |u| ~ 1)
    peclet_number = h * Re
    
    # Recommended time step (most restrictive condition)
    recommended_dt = min(cfl_limit, diffusion_limit)
    
    return (
        h = h,
        cfl_limit = cfl_limit,
        diffusion_limit = diffusion_limit, 
        peclet_number = peclet_number,
        recommended_dt = recommended_dt,
        is_diffusion_limited = diffusion_limit < cfl_limit
    )
end

"""
    find_grid_index(coord::Float64, grid::Vector{Float64}) -> Int

Find the grid index closest to a given coordinate.

# Arguments
- `coord::Float64`: Target coordinate
- `grid::Vector{Float64}`: Grid coordinate array

# Returns
- `idx::Int`: Index of closest grid point
"""
function find_grid_index(coord::Float64, grid::Vector{Float64})
    _, idx = findmin(abs.(grid .- coord))
    return idx
end

"""
    find_special_indices(x::Vector{Float64}, y::Vector{Float64}) -> NamedTuple

Find indices of special locations in the L-shaped domain.

# Returns
Named tuple with:
- `i_x0`: Index where x ≈ 0
- `j_y0`: Index where y ≈ 0  
- `lid_indices`: Indices of moving lid nodes
- `corner_index`: Index of reentrant corner (x=0, y=0)
"""
function find_special_indices(x::Vector{Float64}, y::Vector{Float64})
    # Find zero crossings
    i_x0 = find_grid_index(0.0, x)
    j_y0 = find_grid_index(0.0, y)
    
    # Find moving lid indices (top boundary, x ≤ 0) 
    j_top = length(y)  # Top boundary index
    lid_i_indices = [i for i in 1:length(x) if x[i] <= 0.0]
    lid_indices = [(i, j_top) for i in lid_i_indices]
    
    # Reentrant corner (origin)
    corner_index = (i_x0, j_y0)
    
    return (
        i_x0 = i_x0,
        j_y0 = j_y0,
        lid_indices = lid_indices,
        corner_index = corner_index,
        j_top = j_top
    )
end

"""
    compute_domain_statistics(domain::LShapeDomain{D}) -> NamedTuple where D

Compute comprehensive statistics about the domain and grid.
"""
function compute_domain_statistics(domain::LShapeDomain{D}) where D
    counts = count_domain_nodes(domain)
    
    # Estimate grid spacing (assume uniform grid)
    if D == 2
        nx, ny = size(domain.fluid_mask)
        h_est = 2.0 / (nx - 2)  # Approximate from domain size  
        total_dof = counts.interior
    else # D == 3
        nx, ny, nz = size(domain.fluid_mask)
        h_est = 2.0 / (nx - 2)
        total_dof = counts.interior
    end
    
    # Memory estimates (rough)  
    bytes_per_dof = 8  # Float64
    matrix_memory_mb = (7 * total_dof * bytes_per_dof) / (1024^2)  # 7-point stencil
    vector_memory_mb = (5 * total_dof * bytes_per_dof) / (1024^2)  # Multiple work vectors
    
    return (
        dimension = D,
        grid_spacing = h_est,
        total_dof = total_dof,
        matrix_memory_mb = matrix_memory_mb,
        vector_memory_mb = vector_memory_mb,
        domain_efficiency = counts.interior / counts.total,
        boundary_fraction = counts.boundary / counts.fluid
    )
end

"""
    validate_grid_quality(h::Float64, Re::Float64; warn_threshold::Float64=100.0) -> Bool

Validate grid quality for numerical stability.

# Arguments
- `h::Float64`: Grid spacing
- `Re::Float64`: Reynolds number  
- `warn_threshold::Float64`: Warning threshold for grid Péclet number

# Returns
- `is_acceptable::Bool`: Whether grid quality is acceptable
"""
function validate_grid_quality(h::Float64, Re::Float64; warn_threshold::Float64=100.0)
    metrics = compute_grid_metrics(h, Re)
    is_acceptable = true
    
    # Check grid Péclet number
    if metrics.peclet_number > warn_threshold
        @warn "High grid Péclet number: $(metrics.peclet_number) > $warn_threshold"
        @warn "Consider refining grid or reducing Reynolds number for accuracy"
        is_acceptable = false
    end
    
    # Check if time step will be very restrictive
    if metrics.recommended_dt < 1e-6
        @warn "Very small recommended time step: $(metrics.recommended_dt)"
        @warn "This may lead to slow convergence"
    end
    
    # Check memory requirements (rough estimate)
    if metrics.peclet_number > 1000
        @warn "Very high Péclet number may require special discretization"
        @warn "Consider using upwind schemes or artificial viscosity"
        is_acceptable = false
    end
    
    return is_acceptable
end

"""
    print_grid_info(x::Vector{Float64}, y::Vector{Float64}, Re::Float64)
    print_grid_info(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, Re::Float64)

Print comprehensive grid information for diagnostics.
"""
function print_grid_info(x::Vector{Float64}, y::Vector{Float64}, Re::Float64)
    h = x[2] - x[1]
    metrics = compute_grid_metrics(h, Re)
    
    println("2D Grid Information:")
    println("  Grid size: $(length(x)) × $(length(y))")
    println("  Interior points: $(length(x)-2) × $(length(y)-2) = $((length(x)-2)*(length(y)-2))")
    println("  Grid spacing h: $(@sprintf("%.6f", h))")
    println("  Reynolds number: $Re")
    println("  Grid Péclet number: $(@sprintf("%.2f", metrics.peclet_number))")
    println("  CFL limit: $(@sprintf("%.2e", metrics.cfl_limit))")
    println("  Diffusion limit: $(@sprintf("%.2e", metrics.diffusion_limit))")  
    println("  Recommended Δt: $(@sprintf("%.2e", metrics.recommended_dt))")
    println("  Stability limited by: $(metrics.is_diffusion_limited ? "diffusion" : "CFL")")
    
    validate_grid_quality(h, Re)
end

function print_grid_info(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, Re::Float64)
    h = x[2] - x[1]
    metrics = compute_grid_metrics(h, Re)
    
    println("3D Grid Information:")
    println("  Grid size: $(length(x)) × $(length(y)) × $(length(z))")
    interior_count = (length(x)-2) * (length(y)-2) * (length(z)-2)
    println("  Interior points: $(length(x)-2) × $(length(y)-2) × $(length(z)-2) = $interior_count")
    println("  Grid spacing h: $(@sprintf("%.6f", h))")
    println("  Reynolds number: $Re")
    println("  Grid Péclet number: $(@sprintf("%.2f", metrics.peclet_number))")
    println("  Recommended Δt: $(@sprintf("%.2e", metrics.recommended_dt))")
    
    # Memory estimates for 3D
    dof_estimate = interior_count * 0.6  # Account for L-shape removal
    memory_mb = (7 * dof_estimate * 8) / (1024^2)  # 7-point stencil, Float64
    println("  Estimated memory: $(@sprintf("%.1f", memory_mb)) MB")
    
    validate_grid_quality(h, Re)
end